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Abstract. We study both analytically and numerically the disruptive 
effect of instantaneous gas removal from an embedded cluster. We setup 
a calculation based on the stellar velocity distribution function, to com- 
pute the fraction of stars that remain bound once the cluster has ejected 
J> ■ the gas and is out of equilibrium. We find tracks of bound mass-fraction 

vs star formation efficiency similar to those obtained with N-body cal- 
culations. We use these to argue that embedded clusters must develop 
! high-binding energy cores if they are to survive as bound clusters despite 

q ■ a star formation rate as low as 20% or lower suggested by observations. 
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Q-f 1. Introduction 
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5— ( ' Most if not all stars are formed in embedded clusters and associations, and 

therefore this is likely the predominant mode of star formation contributing to 
the Galactic-field population. To understand how such aggregates form and 
evolve remains a severe challenge however, since no theory exists yet that ac- 
counts in detail for the (presumably) simpler case of the formation of individual 
stars. On the observational front, surveys of star- forming regions suggest that 
the mass-fraction of gas used up at birth does not exceed 10 to 20% of the 
total (Lada 1999). This low star-formation efficiency (sfe) implies that young 
stars or clusters of stars should be embedded in gas. Yet populous clusters with 
ages ^ 1 Myr (e.g. the Orion Nebula Cluster, R136 in 30 Doradus) are already 
void of gas. The fraction of gas left behind after the epoch of star formation 
must therefore be removed quickly (well before the first supernova explodes) 
to reveal a bare cluster. By the time this occurs, most of the cluster must be 
assembled. Hydrodynamic collapse solutions point to the rapid formation of 
stellar cores followed by time-dependent accretion (Boily & Lynden-Bell 1995). 
Stellar masses are then accrued over a period 10 5 — 10 6 years, long before the 
gas-evacuation timescale. Therefore gas-evacuation itself inhibits further star 
formation in the cluster. The gas is driven out by OB stars that culminate the 
star-formation process. Stellar winds from these OB stars yield a momentum 
flux T ~ 8 x 10~ 3 Mq kmsec _1 y r 1 (Churchwell 1999). This is sufficient to blow 
out gas from an 10 3 Mq, 1 pc radius embedded cluster in ~ 10 5 years, which is 
comparable to its crossing time t cr = 2R/a. In addition, the ionising radiation 
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heats the gas to 10 K, causing an overpressure and expansion at the sound 
velocity « 10 km/s. We want to establish how significant mass-loss affects the 
structure of a cluster, in particular what fraction of the stars remain bound. 

This contribution marks the start of a research programme aiming at iden- 
tifying and quantifying the fundamental physical processes responsible for the 
formation of open clusters. This programme is inspired by the finding from high- 
precision iV-body computations that open clusters can form readily despite low 
star-formation efficiencies, in contradiction to previous results (Kroupa, Aarseth 
& Hurley 2001). Three hypothesis are raised that may explain these A-body 
results, (i) Either two-body encounters during the expansion phase of the clus- 
ter after gas expulsion condense a part of the radial flow to a bound entity, or 
(ii) encounters between the ubiquitous binary systems during the radial outflow 
affect the condensation, or (iii) two-body and binary-binary encounters before 
gas-expulsion (i.e. during the embedded phase) cause the segregation of a tightly 
bound core that resists expansion when the gas leaves the system. 

In order to weigh the relative contributions from each of points (i)-(iii), 
all drawn from particle-particle (collisional) evolution, we seek first to isolate 
the salient features of collisionless, smooth open cluster evolution. We take an 
analytic approach to the problem, based on the velocity distribution function of 
stars. By applying a new, fast and self-consistent iterative method for computing 
the fraction of bound stars, we show that to account for observations of low 
sfe, clusters must develop strongly bound cores to avoid dissolving completely 
owing to mass loss. In addition to our analytical approach we use numerical 
calculations to assess its range of applicability. 

2. The problem 

A cluster of stars forms converting a fraction e of the gas into stars in the 
process. Stellar winds or a supernova event blow/s out the remainder on a 
timescale r <C t cr . What fraction of the stars remain to form a bound cluster ? 
How is the equilibrium profile linked to the initial mass density ? Since the star 
formation epoch will have lasted as long as or longer than a cluster dynamical 
time t CI , the system as a whole is close to virial equilibrium before gas is expelled. 

Hills (1980) argued that equilibrium, self-gravitating stellar systems would 
expand or even dissolve if half the mass is lost very quickly: stars then preserve 
their kinetic energy established under a deeper potential well, hence may escape 
if their binding energy becomes positive. We write the total energy 



where M, R are the mass and radius of a uniform-density spherical distribution 
of mean square velocity (v 2 ), and k = 3/5 for the particular geometry and 
density profile considered. Before gas-expulsion the total mass in gas and stars 
is M = M init = M gas + M*, while after the gas is expelled M = M± but the 
velocity dispersion, (v 2 ) = KGM [nit /R, from the scalar virial theorem. Thus we 
find a solution for E > if M+ < Mj n i t /2, i.e. as the mass is reduced by 50% or 
more (e < 1/2), the remaining system has zero or positive energy globally. 
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But if the law of averages applies, nothing stops a hard-bound core forming 
at the expense of an expanding, loose envelope of stars, while E > for the 
system as a whole. We need to find out under which circumstances this will 
occur. To do this we introduce an iterative procedure based on the stellar 
distribution function in order to determine the end-product fraction of bound 
stars. 



3. Distribution function: iterative scheme 



As a simplifying assumption we take the sfe to be independent of position in the 
cluster. Hence e = constant (see Adams 2000 for a different approach). Thus 
gas and stars are initially mixed in the same proportion throughout, however 
note that e does not fix the profiling of stellar masses with radius, rather the 
total mass of gas made into stars : 

V: (number of stars in mass bin i) x (stellar mass m,) 

=^ = constant 

^gas 

at any radius r. Both gas and stars are distributed according to the same 
distribution function which for a spherical system depends only on energy, f(E). 
The total system mass is then 

M = M* + M gas = e" 1 M* = \ f(E)dE, 

Je c 

where \E C \ is the maximum binding energy and with a suitable normalisation of 
f(E). The cluster potential 4>(r) may be written 

rr A lt ru OM 

cP(r;M,R)= — ^Gp(w)w 2 dw = y(r/R) -—, (1) 
Joo u z Jo R 

where y(x) is a dimensionless function. Keeping R constant while removing the 
gas instantly so the potential includes stars only, we have 



0(r;M,12)-& = e6 (2) 
and hence the fraction of stars at radius r which now have positive energy is 



f(v,r/R)v dv 

I - /. ^ — , (3) 

f(v,r/R)v dv 

with v e the (local) escape velocity computed for <f> (before gas-expulsion), and 
similarly for v e il obtained for 0*. 

Thus for a given d.f. and e, we may compute the quantity 1 — f e at all radii. 
Note that should the sfe be a function of radius, the simple renormalisation (2) 
leading to 0* would not apply, however (||) may still be computed if the stars' 
potential is given and v e known from the initial gas + stars mixture. The key 
step is to adjust the potential 0* itself, since potential and velocity field do 
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not match any longer. This is normally computed using N-body integration to 
take into account the full star-star interactions, and the redistribution of energy 
between them in the time-dependent potential. If we thought that stars acquire 
only little energy during the time that they escape, then the fraction of stars 
remaining might be computed as follows. 

Since the fraction (||) of positive-energy stars is known at all r, we recom- 
pute the gravitational potential counting only stars with E < at each radius. 
Neglecting dynamical evolution, the cluster radius R is unchanged and hence 
the potential can be recomputed by integration, from oo, inwards. Once the 
new potential is known, a fraction of the remaining stars will again be unbound 
by virtue of the stars lost during the previous iteration. We therefore repeat the 
procedure until finally the cluster mass converges to a finite quantity, in which 
case no more stars escape and the original distribution function is depleted from 
all escapers in a self-consistent manner. 

We consider three cases in detail, then turn briefly to numerical calculations. 

3.1. Power- law d.f. 

The case f(v) oc v@~ 2 , where (5 is constant, provides a useful illustrative starting 
point. The d.f. is truncated at the local escape velocity. Then the velocity dis- 
persion <t 2 oc v 2 = 2(f)(r) at each radius. We find, on substituting p(w) — > ep(w), 
the new potential given by @ and i> e ,*(r) = ^/2</>*(r) = e 1//2 f e (r). Equation (||) 
becomes (with u = v/v e ) 



f(v,r/R)v 2 dv f f( u )u 2 du 
1 _ U = = = i _ e 03+D/2 > (4) 

/ f(v,r/R)v 2 dv f f( u )u 2 du 
J o Jo 

where we took f3 > —1. Thus f e = e^ +1 ^ 2 is independent of radius. Repeating 
the procedure to take account of the positive-energy stars, we substitute e — > 
f e • e, etc, so that after n iterations the net mass of bound stars, M*, becomes 



M b = e [w+mr . e p+i)/2]»-\ .. e (/3+i)/2 eM = n n =i ( e p+i)/2]^ M ^ (5) 

(and similarly for the stellar density p*[r] wrt p[r\). The multiplicative operator 
II in (|5|) leads to a non-zero (positive) value as n — > oo only when /3 < 1, since 
e < 1. All d.f.'s with (5 > 1 lead to cluster disruption, because the high- velocity 
range of the d.f. is too densely populated, leading to catastrophic stellar loss 
after the expulsion of any amount of gas. 

3.2. Plummer model 

We wish to compare our basic result (|5|) to a standard fit to globular clusters. 
We consider the Plummer model, where f(E) oc (—E) 7 ^ 2 . For this case the total 
system mass is finite but infinite in extent; the velocity dispersion maximises at 
the centre, as the density. The velocity d.f., f(v) oc (1 — (v/v e ) 2 ) 7 ^ 2 (see Spitzer 
1987); writing p(e) = 105 - 1210e + 2104e 2 - 1488e 3 + 384e 4 , we find 
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(1 — U 2 ) 7 / 2 U 2 du I/O R / \ inr ■ -1 1/2 

i- /e = ^ — ! = i + e/vr ^ p( l; 1058111 e/ . (6) 



(l- U 2 ) 7 / 2 u 2 du 

o 

We may repeat the procedure until 1 — f e — > and no additional stars are lost. 
We were not able to express the resulting expression in simple form, however we 
note that, as in the first case, the solution is independent of radius, hence the 
density profile is simply renormalised at each radius. Therefore the potential 4>* 
and escape velocity follow from (||). 

3.3. Hernquist profile 

To contrast with the smooth-density Plummer solution, we consider a peaked 
Hernquist (1990) profile. The density p and one-dimensional velocity dispersion 
cr vary radially according to 

, Mil Mil 1 4>(x) 

p(r or x — — 



a (x) = 4>(x) x (1 + x) < In 



2tt r (r + r c ) 3 27rr 3 x (x + l) 3 2irGr^ x(x + l) 2 ' 
l + x 1/4 1/3 1/2 1 



(l+x) 4 (l + x) 3 (l + x) 2 l + x. 

with r c a free length fixing the point of the power-law turnover, and x = r/r c . 
The velocity d.f., f(v), is only constrained locally by <r(r); we thus have the 
freedom to choose any profile satisfying a(r). We set f(v) oc v 2 exp(— v 2 /2a 2 ), a 
Maxwellian profile. This is found to give stable equilibria in N-body calculations 
(Hernquist 1993). Inserting this in (^) yields 

V^e"* - erf(V2^) 

where the dimensionless potential \& = 4>/cr 2 and erf(x) is the error function. 
Note that although the dispersion a 2 oc <p as before, here the fraction of positive- 
energy stars depends on the local potential, and hence it is a function of radius. 
To compute the net fraction of bound star we must therefore recompute the 
potential numerically for each evaluation of f e in (|7j). This poses no problem 
since (f>(r — ► do) — > 0, and the density is known at each step (though it is no 
longer a Hernquist profile). 



4. Results for three illustrative cases 

Our results are shown in Fig. 1 for the three cases discussed above. The solid 
lines show the solutions for the power-law d.f. with (3 = —3/4 and 0. Note that 
in either case the fraction M^/M* of stars that remain bound is not dropping 
to zero until e itself is zero. The case (3 = corresponds to a flat distribution, 
f(v)v 2 dv = constant xdv. Assuming only that no star has velocity greater than 
the local escape velocity, the solution (||) decreases only marginally faster than 
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e as e —* 0. By contrast, a Plummer or Maxwellian-Hernquist model shows a 
sudden drop to a null fraction of bound stars for finite sfe : for Plummer models 
we obtain a critical value e ~ 0.445, while for the Hernquist model the fraction of 
bound stars exceeds 5% or so until e ~ 0.28. In both cases the iterative scheme 
converged to machine accuracy. 

Note that for the two self-consistent d.f.'s discussed here, the power-law 
solutions provide an illustrative description of the survival rate as e approaches 
a critical value. The more robust Hernquist model favours low- velocity stars 
(near the centre a — > 0) and hence is better fitted with the solution (3 = —3/4 
in the range 0.5 < e < 1. Figure 2 shows for this case that the fraction of 
bound stars f e remains larger near the core (see Fig. 2b). As a result, the 
initial density profile becomes steeper with radius. Since we have only made a 
selection by energy, the expectation is that the bound system indeed should be 
more peaked than initially. 

.1 I I I | I I I I | I I I I | I I I I | I I I I | I I I I | I I I I | I I I I | I I I I | I I I I | I I I L 



1 - 




-■ i ■ ■ I ■ ■ ■ i I ■ ■ ■ ■ I ■ ■ ■ ■ I i ■ ■ ■ I ■ ■ ■ ■ I ■ ■ i ■ I ■ ■ ■ ■ I ■ ■ i ■ I ■ ■ ■ ■ I ■ i ■ r 



0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 
star formation efficiency 

So far we only considered a mapping of static configurations under the se- 
lection of particles according to (||). In reality stars must orbit in space before 
leaving the system, and they may exchange gravitating energy with the back- 
ground system in the process. N-body calculations are of some help here. N- 
body studies with time-dependent potential (Lada, Margulis &: Deardorn 1984) 
or a star-gas mixture (Geyer & Burkert 2000) also find a critical sfe below which 
clusters dissolve. We decided to conduct our own N-body calculations with a 
collisionless grid code (Fellhauer et al. 2000) and 100,000 particle equilibrium 
Plummer models. Our analytic approach gives a critical value for survival of 
e = 0.4448. We therefore setup two N-body calculations with sfe = 50% and 
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40%. The fraction of stars which remain bound in each case brackets the results 
derived from (g) (see Fig. 1, large open circles). Notably we find no indication 
that any stars remain for the case where e = 0.40, in agreement with the findings 
by others (Lada, Margulis & Deardorn 1984; Geyer & Burkert 2000) for similar 
setup 
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Figure 2. Initial and final density profiles for (a) Plummer and (b) 
Hernquist models for two values of the sfe e. The run of f e , i.e. the ratio 
of density to the initial density, is also shown as a function of radius. 
Note how the Hernquist profile becomes steeper for small e. 



5. Time-evolution and other considerations 

To progress further, we note that only two timescales are important to the 
problem, namely the cluster crossing time, i cr , and the gas removal timescale, 
r. We may understand the dynamics by comparing r to t CT . 

5.1. Collisional evolution, but short r 

The key lies with the redistribution of kinetic energy between the stars. This 
can be achieved on a short timescale by direct collisions or close encounters, 

when the Safronov cross section d 2 + N y 3 ^ is large (here d 2 is a star's ge- 
ometric cross section and N the number of stars). This is especially true for 
small-N open clusters with binaries and multiple stars, in which case d is set 
equal to their semi-major axis. In this situation collisional effects are never neg- 
ligible and hence when r ~ t CT or longer, the situation is not one of equilibrium, 
and evolution must be tackled numerically. Kroupa, Aarseth & Hurley (2001) 
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evolved an embedded Plummer model with a high-precision direct-summation 
iV-body code and delayed, but near-to instantaneous, gas-removal. They find 
a fraction ~ 30% of stars remain despite a low sfe of e ~ 0.3. The results are 
similar to the analytical results for the Hernquist model (Fig. 1, black square), 
while we would have expected complete dissolution for collisionless evolution of 
a Plummer model. We are lead to the conclusion that it is the compact core that 
develops as a result of two-body relaxation during the embedded phase before gas 
expulsion that leads to more robust clusters (see also Section 1). 

5.2. Collisionless evolution, but long r 

Large-N systems possess a long two-body relaxation time t co i ~ ToWjN ^ cr " When 
t C r *C t <C i C ob two-body effects may be neglected. In this case the age of the 
cluster may not have allowed for global two-body relaxation, yet significant mass 
removal will have occurred over several stellar orbits concentrated around the 
centre and these orbits evolve adiabatically. Since adiabatic evolution is a re- 
mapping of an orbit to itself, no stars on such orbits are lost. Thus for finite 
or large r, we anticipate the survival rate of clusters to be intimately linked to 
their properties at birth, such as what family of orbits are present initially. 

The future of this programme will see additional high-accuracy iV-body 
computations being performed to address how relatively important the three 
hypothesis raised in Section 1 are for the formation of bound clusters. Specif- 
ically, the formation of sub-condensations as a result of encounters during the 
radially expanding flow will be addressed in detail. 
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